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Transitions between "glassy" local minima of a model free- 
energy functional for a dense hard-sphere system are studied 
numerically using a "microcanonical" Monte Carlo method 
that enables us to obtain the transition probability as a func- 
tion of the free energy and the Monte Carlo "time". The 
growth of the height of the effective free energy barrier with 
density is found to be consistent with a Vogel-Fulcher law. 
The dependence of the transition probability on time indi- 
cates that this growth is primarily due to entropic effects 
arising from the difficulty of finding low-free-energy saddle 
points connecting glassy minima. 
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The dynamic behavior of supercooled liquids near the 
glass transition Q is one of the most enigmatic problems 
of condensed matter physics. The most dramatic fea- 
ture of the dynamics near the glass transition in so-called 
"fragile" systems Q is an extremely rapid growth of the 
relaxation time r which is reasonably well-described by 
the Vogel-Fulcher law j|, r oc exp[C/(T - T )], where 
To < T g , the conventionally defined glass transition tem- 
perature at which the viscosity attains a value of 10 13 P. 
The apparent divergence of r at To has led to specula- 
tions about the possibility of a true thermodynamic tran- 
sition at this temperature. This is also suggested by the 
observation Q that the temperature Tjc (the so-called 
Kauzmann temperature) at which the entropy difference 
between the supercooled liquid and the equilibrium crys- 
talline solid extrapolates to zero is very close to To. The 
closeness of Tq and Tjc suggests that the growth of the 
relaxation time near the glass transition is primarily en- 
tropic in origin. Heuristic arguments that attempt to 
relate the Vogel-Fulcher law to entropic effects have been 
proposed by several authors J^||. However, we are not 
aware of any calculation that provides an explicit demon- 
stration of such effects in simple model liquids. 

In this Letter, we describe the results of a numerical in- 
vestigation that provides direct evidence for an entropic 
origin of the growth of the relaxation time in a dense hard 
sphere fluid near the glass transition. Our computations 
are based on a model free energy functional for the 
hard sphere system. We use a novel "microcanonical" 
Monte Carlo (MC) method to study transitions between 
different "glassy" minima of a discretized version of this 



free energy functional. We determine the probability of 
transition from one minimum to another as a function of 
the free energy increment A/ (the excess free energy per 
particle measured from that at the original minimum) 
and MC "time" t. This allows us to define an effective 
barrier height that depends weakly on t. We find that 
the growth of this effective barrier height with increasing 
density is consistent with a Vogel-Fulcher form appro- 
priate for a hard-sphere system. Our numerical results 
about how the dependence of the effective barrier height 
on t changes as the density is increased indicate clearly 
that the growth of the barrier height (and the conse- 
quent growth of the relaxation time) is primarily due to 
entropic effects arising from an increase in the difficulty 
of finding low free-energy paths ("saddle points") that 
connect one glassy local minimum with another. 

The free energy functional used in our study is of the 
form proposed by Ramakrishnan and Yussouff @ : 



F[p]=F l \p Q ]+k B T 



dr{p(r) ln(p(r)/p ) - Sp(r)} 



(1/2) / dr / dr'C(\r-r'\)Sp(r)8p(r') 
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where Sp(r) = p(r) — po is the deviation of the time- 
averaged local number density p(r) from its value po in 
the uniform liquid state, Fi is the free energy of the uni- 
form liquid, T is the temperature, and C(r) is the di- 
rect pair correlation function Q of the uniform liquid at 
density pq. C(r) is expressed in terms of the dimension- 
less density n* = poa 3 (ex is the hard sphere diameter) 
through the Percus-Yevick approximation Q which is ex- 
pected to be adequate if po is not very high. 

The discretized version of this free energy functional 
exhibits |9| a large number of "glassy" minima (local 
minima of F at which the density is inhomogeneous but 
aperiodic) for n* > n*j where rfj ~ 0.85 is the density 
at which equilibrium crystallization occurs. Numerical 
studies [ |lO|]ll| ] of Langevin equations appropriate for this 
system show that the dynamic behavior is governed by 
thermally activated transitions among these glassy min- 
ima if n* exceeds a "crossover" value that is close to 
0.96. The time scales for such transitions were estimated 



from a standard MC method in Ref. |12| and found to 
rapidly increase with increasing density. Here, we have 
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used a different numerical method that is more efficient 
than the canonical MC method and provides information 
about the origin of the growth of the time scale of these 
thermally activated transitions. The increase of this time 
scale may be due to one or both of two factors: (a) an 
increase of the height of the lowest free-energy barrier 
that separates a glassy minimum from the others; and 
(b) an increase of the difficulty in finding low-free-energy 
paths to other minima. Considering the free energy func- 
tional as an effective Hamiltonian for the system, these 
two factors may be called "energetic" and "entropic" Q , 
respectively. The canonical MC method does not provide 
much information about the relative importance of these 
two factors in the observed growth of the time scales. 
As described below, the numerical method used in this 
study allows us to distinguish between energetic and en- 
tropic effects. It also allows us to follow the growth of the 
barrier-crossing time scale over about ten decades, which 
would not be possible in a canonical MC calculation. 

We discretize our system on a cubic lattice of size L 3 
and mesh constant h with dimensionless density variables 
defined as p L = p(r.i)h 3 . Periodic boundary conditions 
are used and the constraint that the sum of the variables 
Pi must be a constant N, the number of particles in the 
sample, is enforced during the simulation. We define a 
dimensionless free energy per particle, f[p], as: 

f[p] = (3F[p]/(n*L 3 a 3 ) (2) 

where a is the ratio h/a, and (3 — 1/fcsT. 

Our numerical method, which may be called "micro- 
canonical" MC if the free energy functional is consid- 
ered to be an effective Hamiltonian, involves the following 
steps. Each run is started from a glassy local minimum of 
the free energy. We choose a trial value of what we call 
the free energy increment Af and then perform a MC 
simulation in which we sweep the sites i of the lattice se- 
quentially. At each step and site, we pick another site j 
at random from the ones that lie within a distance a from 
the site i. We then attempt to change the values of pi and 
Pj to p(pi + pj) and (1 — p)(pi + pj), where p is a random 
number distributed uniformly in [0,1]. The atempted 
change is accepted only if the dimensionless free energy 
after the change is less than F m i n + NAf where F m i„ is 
the dimensionless free energy [3F at the minimum where 
the simulation is started. This procedure generates a 
random sampling of configurations whose free energy lies 
within NAf of that of the glassy minimum under con- 
sideration. The simulation proceeds up to a maximum 
"time" t m , of MC steps per site. At regular intervals 
along the evolution of the system, we use a minimiza- 
tion procedure || to determine whether the system has 
moved to the basin of attraction of a different local min- 
imum of the free energy. Obviously, if NAf is smaller 
than the lowest free-energy barrier between the starting 
minimum and any other minimum, the system remains 



in the basin of attraction of the starting minimum. As 
Af is increased, one begins to find transitions to other 
accessible minima, that is, minima that the system can 
find within a time t < t m , which are separated from the 
initial one by a barrier of height less than NAf. Repeat- 
ing this procedure a number of times (typically 10-20) 
for a fixed set of values of n*, Af and t m , we obtain 
P(n*, Af, t), the probability of a transition to a different 
minimum within time t for free energy increment Af, 
as the fraction of the number of runs in which a transi- 
tion is found. This probability is calculated for a suitable 
range of values of n* , Af, and t, and the whole procedure 
is repeated for several glassy minima of the free energy 
(see below). We define a "critical" value, Af c (n*,t), of 
the free energy increment as the value of Af for which 
P(n*, Af, t) — 0.5. Clearly, NAf c represents an effective 
barrier height for transitions to other local minima. This 
is the quantity that we use to present our results. 

We have used two sizes, L = 15 and L = 12. In the first 
case we have taken a = 1/4.6 so that L and a are incom- 
mensurate with a close-packed lattice and no crystalline 
minimum of the free energy is found. The total number of 
inhomogeneous minima is then about 10 and all of them 
exhibit glassy structure as determined by the two-point 
correlation function of the local density. The minima we 
have used as our starting point in this case were also 
used in Ref. |12|. These are the minima to which the 
system moves TO] in the course of its time evolution un- 
der Langevin dynamics |l0| when it is started from the 
uniform liquid state. For L — 12 we took a = 0.25 so 
that the sample is commensurate. It admits a crystalline 
minimum that has the lowest free energy for the values of 
n* considered here. The number of glassy minima is sub- 
stantially larger (about 30) in this sample. Out of those 
we chose a few with structure similar to that of the min- 
ima of the L = 15 sample. For both cases, the minima 
found at lower densities were "followed" to higher densi- 
ties by running the minimization program at the higher 
density using the lower density configuration (which is 
of course not a minimum at the higher density) as the 
starting point. The values of t m are 15,000 for L = 15 
and 8,000 for L = 12. The transition probability was 
calculated at time intervals of 5,000 in the first case, and 
2,000 in the second case. In both cases, the density range 
covered was 0.94 < n* < 1.06. Higher values of n* were 
not considered because the Percus-Yevick approximation 
then becomes || inaccurate. 

Typical results for P are shown in Fig.l where data 
for L = 12 and n* — 1.04 are plotted for four different 
values of t. The value of Af was incremented in steps 
of 0.05, which is also the estimated uncertainty in the 
determination of Af c . The transition probability grows 
from zero as Af is increased, and eventually saturates 
at one for sufficiently large values of Af. For a fixed 
value of Af, the transition probability increases as t is 
increased: transitions to other minima are more likely 
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if the system is allowed to explore a larger number of 
configurations. Since P is an increasing function of both 
Af and t, A/ C (n*,i) (the value of A/ where P=0.5, as 
defined above) decreases as t is increased. In agreement 
with the previously observed Jl2[ growth of the barrier- 
crossing time scale with n*, we find that A/ c is an in- 
creasing function of n* . 

The conclusion that entropic effects play a crucial role 
in the growth of the effective height of the free energy bar- 
riers stems from the observation that the i-dependence 
of A/ c becomes stronger as n* is increased (see Fig. 2). 
It is intuitively obvious that the i-dependence of A/ c is 
closely related to the probability of finding a path ( "sad- 
dle point" ) that connects the starting minimum to a dif- 
ferent one. If such paths were relatively easy to find, 
then the transition probability would be insensitive to the 
value of i as long as it is not very short. If, however, paths 
to other minima are few, a large number of configurations 
have to be explored before one of them is found. The t- 
dependence of A/ c would then be more pronounced and 
extend to larger values of t. To make the idea more con- 
crete, we ignore the correlations (which are short-range 
in time) among the configurations generated in a MC run 
and assume that they represent t independent samplings 
of configurations with free energy less than F m i n + NAf. 
Let us also assume that the system does not return to 
the basin of attraction of the starting minimum after a 
transition to a different basin of attraction. We find that 
a return to the original basin of attraction is indeed very 
rare. The transition probability may then be estimated 
as P(n*,Af,t) = 1 - [l-p(n*,A/)]* ~ 1 - exp(-fp), 
where p{n*,Af) < 1 is the probability that a ran- 
domly chosen configuration with j3F < F m i n + NAf 
belongs in the basin of attraction of a different mini- 
mum. One expects p to be zero if Af < A/o(n*) where 
NAfo is the height of the lowest free energy barrier, and 
p = g{n*, Af — A/o) for Af > A/o where g(n*, x) grows 
continuously from zero as x is increased from zero. Com- 
bining this with the definition of Af c , we obtain the re- 
lation g(n*, A/ C (n*, t) - Af Q (n*)) = ln2/t. Our observa- 
tion that the difference Af c (n*, t\) — Af c (rb*,t%) for fixed 
t\ < <2 increases with n* then leads to the conclusion that 
the function g(n*,x) decreases (i.e. the difficulty of find- 
ing paths to other minima increases) as n* is increased 
at fixed x. 

The observed ^-dependence of A/ c for all values of n* 
and all the minima in our study is well-represented by 



Af c (n*,t)=Af Q (n*) + c(n*)t- 



(3) 



with a in the range 0.25 — 0.40. Typical fits to this form 
with a = 0.35 for two minima with L = 15 and L = 12 
are shown in Fig. 2. The values of A/o obtained from such 
fits with a fixed value of a are nearly independent of n* , 
but exhibit a dependence on the value of a, varying be- 
tween and 0.5 for the L = 15 minimum, and between 1.3 



and 1.5 for the L — 12 minimum of Fig. 2. The quantity 
c(n*) increases with n*. These results correspond to the 
function g(n* : x) having the form g(n* : x) ~ A(n*)x 1 / a 
with A(n*) decreasing with increasing n* . We conclude 
from these observations that the growth of the effective 
barrier height with increasing n* is primarily due to an 
entropic mechanism associated with an increase of the 
difficulty in finding low-lying saddle points that connect 
different glassy local minimum of the free energy. This 
conclusion is consistent with the canonical MC results 
of Ref. [jl2j where we found that while the time scale of 
transitions between minima increases dramatically with 
n* , the free energy increment at the transition point re- 
mains essentially independent of n*. This implies that 
the probability of finding a path to other minima for a 
fixed value of the free-energy increment decreases as n* 
is increased. 

Our results for the dependence of A/ c on n* are con- 
sistent with the Vogel-Fulcher law || which assumes the 
following form |14| for our system: 



A/ C (n*) = a + b/(n* c - n*), 



(4) 



where a and b are constants and n* is expected to be 
close to the random close packing density n* cp ~ 1.23. 
There is some ambiguity in trying to fit our data to this 
form because our values of A/ c depend weakly on the 
time t. However, the value of n* obtained from fits of 
our data for Af c (n*,t) to Eq.(||) with fixed a is nearly 
independent of t. This is consistent with the form of 
Eq.(|) if a = A/o , b cx t~ a , and c cx l/(n* - n*). A/ is 
indeed nearly independent of n*, and the ^-dependence 
of b and the n*-dependence of c are in agreement with 
the other two conditions. For the L = 15 case, we can 
fit the data for A/ c at t = 15,000 to the form of Eq.(f§) 
with a = (A/o = 0). The best fit, shown in Fig. 3, 
corresponds to n* c — 1.225, very close to the expected 
result. The best fit to the L = 12 data with a ~ 1.0 
(the difference between the values of A/o for the L = 12 
and L = 15 minima is about 1.0) also yields a similar 
value of n*. So, we conclude that the observed growth of 
the effective barrier height is consistent with the Vogel- 
Fulcher form. The increase in the effective barrier height 
as n* is increased from 0.94 to 1.06 is about 25fcsT, cor- 
responding to a growth of the characteristic time scale of 
about ten orders of magnitude. Thus, the range of time 
scales covered in our study is comparable to that used in 
Vogel-Fulcher fits of experimental data, and much wider 
than what can be achieved in standard MC or molecular 
dynamics simulations. 

In summary, our study shows that the density- 
dependence of the characteristic time scale for transi- 
tions between glassy local minima of the free energy of a 
dense hard sphere system arises primarily from entropic 
effects associated with the difficulty of finding low-lying 
paths that connect such minima. The observed growth 
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of the time scale is quantitatively consistent with the 
Vogel-Fulcher form. To our knowledge, this is the first 
explicit demonstration of entropic effects leading to a 
Vogel-Fulcher-type growth of relaxation times in a simple 
model liquid. The same behavior should occur in other 
simple liquids where C(r) is similar Q. 
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FIG. 3. The dependence of Af c (n*,t = 15000) on n* for L 
— 15. The dashed line shows the best fit to a Vogel-Fulcher 
form (see text). 
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FIG. 1. The transition probability P (see text) as a func- 
tion of the free-energy increment A/ for four values of the 
time t. The data shown are for a L — 12 minimum at n* — 
1.04. The values of A/ c are indicated by the filled circles. 



FIG. 2. Plots of A/c, obtained for a L = 15 minimum, 
against (t/1000) -0 ' 35 for three values of n*. The dashed lines 
are the best straight-line fits. Similar plots for a L = 12 
minimum are shown in the inset. 
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